This document was rendered at 2019-06-20 23:44:17
Figure 1. Cartoon of graph we will build manually via its adjacency matrix.
adj_mat = rbind(c(0,3,1,0,0,0,0,0,0), #A's neighbors
c(3,0,2,2,0,0,0,0,0), #B's neighbors
c(1,2,0,0,2,0,0,0,0), #C's neighbors
c(0,2,0,0,1,0,1,1,0), #D's neighbors
c(0,0,2,1,0,2,0,2,0), #E's neighbors
c(0,0,0,0,2,0,0,0,0), #F's neighbors
c(0,0,0,1,0,0,0,1,0), #G's neighbors
c(0,0,0,1,2,0,1,0,1), #H's neighbors
c(0,0,0,0,0,0,0,1,0) #I's neighbors
)
rownames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
colnames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
# Convert adjacency matrices to igrpah objects for all three graphs.
ig = graph.adjacency(adj_mat, mode="undirected", weighted=TRUE, add.colnames = "name")
print(ig)
## IGRAPH 419525a UNW- 9 12 --
## + attr: name (v/c), weight (e/n)
## + edges from 419525a (vertex names):
## [1] A--B A--C B--C B--D C--E D--E D--G D--H E--F E--H G--H H--I
# Load the Miller2015_Heparin dataset
data(Miller2015_Heparin)
data(heparin_metabolites)
data = Miller2015_Heparin
# Only include metabolites that are present in >90% reference samples.
fil.rate = anno$`Times identifed in all 200 samples`/200
names(fil.rate) = rownames(anno)
data = data[which(fil.rate>0.90), ]
#dim(data)
# Remove any metabolites where any profile has a z-score > 1000. These are likely imputed raw values that were not z-scored.
rmMets = names(which(apply(data, 1, function(i) any(i>20))))
if (length(rmMets)>0) {
data = data[-which(rownames(data) %in% rmMets),]
}
#dim(data)
# Get the diagnoses associated with the samples
data(Miller2015_diags)
# Get data from all patients with Argininemia
arg_data = data[,which(diagnoses$diagnosis=="Argininemia")]
# Add surrogate disease and surrogate reference profiles based on 1 standard deviation around profiles from real patients to improve rank of matrix when learning Gaussian Markov Random Field network on data.
arg_data = data.surrogateProfiles(arg_data, 1, FALSE, TRUE, ref_data = data[,which(diagnoses$diagnosis=="No biochemical genetic diagnosis")])
#dim(arg_data)
# Learn a Gaussian Markov Random Field model using the Graphical LASSO in the R package "huge".
# Select the regularization parameter based on the "STARS" stability estimate.
# require(huge)
#This will take 30 seconds - 1 minute.
arg = huge(t(arg_data), method="glasso")
plot(arg)
# This will take 3-5-ish minutes.
arg.select = huge.select(arg, criterion="stars")
plot(arg.select)
# This is the regularization parameter the STARS method selected.
print(arg.select$opt.lambda)
# This is the corresponding inverse of the covariance matrix that corresponds to the selected regularization level.
arg_icov = as.matrix(arg.select$opt.icov)
# Remove all "self" edges, as we are not interested in self-relationships.
diag(arg_icov) = 0
rownames(arg_icov) = rownames(arg_data)
colnames(arg_icov) = rownames(arg_data)
# Convert adjacency matrices to an igraph object.
ig_arg = graph.adjacency(arg_icov, mode="undirected", weighted=TRUE, add.colnames = "name")
print(ig_arg)
Run the following code, then go to the directory, and open all diffusionEventMovie*.png files all at once. Starting from the first image, view how the probability diffusion algorithm works to diffuse 100% probability to the rest of the graph. Be sure to pay attention to the recursion level listed in the title of each image, to imagine where in the call stack the algorithm is at the captured time the image was generated.
# Set some global parameters for the Probability Diffusion Algorithm.
p0=0.0
p1=1.0
thresholdDiff=0.01
G=vector(mode="list", length=length(V(ig)$name))
G[1:length(G)] = 0
names(G) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
startNode = "A"
visitedNodes = startNode
# Probability diffusion truncates at
thresholdDiff=0.01
coords = layout.fruchterman.reingold(ig)
V(ig)$x = coords[,1]
V(ig)$y = coords[,2]
# Global variable imgNum
imgNum=1
G_new = graph.diffuseP1Movie(p1, startNode, G, visitedNodes, ig, recursion_level=1, output_dir = "/Users/lillian.rosa/Downloads/CTD/vignette/diffusion_event")
# Inherited probabilities across all nodes should add to 1.
sum(unlist(G_new))
## [1] 1
# Which node inherited the highest probability from startNode?
G_new[which.max(G_new)]
## $B
## [1] 0.40625
Now, delete all diffusionEventMovie*.png files from your current directory, and run the following code. View the new image stack in the same way we did previously.
# Now let's see how the probability diffusion algorithm diffuses probability after B has been "stepped" into.
visitedNodes = c("A", "B")
startNode = "B"
imgNum=1
G_new = graph.diffuseP1Movie(p1, startNode, G, visitedNodes, ig, 1, output_dir = "/Users/lillian.rosa/Downloads/CTD/vignette/diffusion_event2")
# Inherited probabilities across all nodes should add to 1.
sum(unlist(G_new))
## [1] 1
# Which node inherited the highest probability from startNode?
G_new[which.max(G_new)]
## $D
## [1] 0.2791667
Sometimes the startNode is “stranded” by a bunch of visited nodes. The diffusion algorithm diffuses “through” visited nodes, so that nodes in the same connected component can be prioritized over nodes in a different connected component, or “island nodes” (e.g. in the below code snippet, “I” is an island node). This only works currently for nodes 2 hops away from the current startNode, however.
adj_mat = rbind(c(0,1,2,0,0,0,0,0,0), #A's neighbors
c(1,0,3,0,0,0,0,0,0), #B's neighbors
c(2,3,0,0,1,0,0,0,0), #C's neighbors
c(0,0,0,0,0,0,1,1,0), #D's neighbors
c(0,0,1,0,0,1,0,0,0), #E's neighbors
c(0,0,0,0,1,0,0,0,0), #F's neighbors
c(0,0,0,1,0,0,0,1,0), #G's neighbors
c(0,0,0,1,0,0,1,0,0), #H's neighbors
c(0,0,0,0,0,0,0,0,0) #I's neighbors
)
rownames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
colnames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
# Convert adjacency matrices to igrpah objects for all three graphs.
ig = graph.adjacency(adj_mat, mode="undirected", weighted=TRUE, add.colnames = "name")
print(ig)
## IGRAPH f620114 UNW- 9 8 --
## + attr: name (v/c), weight (e/n)
## + edges from f620114 (vertex names):
## [1] A--B A--C B--C C--E D--G D--H E--F G--H
adjacency_matrix = list(adj_mat)
# Now let's see how the probability diffusion algorithm diffuses probability after B has been "stepped" into "C" and then "A". As you can see, startNode "A" is surrounded by visited nodes "B" and "C". It needs to be smart enough to weigh "E" and "F" before "D", "H", "G" and "I".
visitedNodes = c("B", "C", "A")
startNode = "A"
G_new = graph.diffuseP1(1.0, startNode, G, visitedNodes, 1, verbose=TRUE)
## [1] " prob. to diffuse:1.000000 startNode: A, visitedNodes: B, C, A"
## [1] " Weighted_edges.sum=1.000000"
## [1] " added diffused probability 1.000000 to child #1: E"
## [1] " subtracted 0.500000 from startNode's neighbor #1: E and sent to its own neighbors."
# Inherited probabilities across all nodes should add to 1.
sum(unlist(G_new))
## [1] 1
# Which node inherited the highest probability from startNode?
G_new[which.max(G_new)]
## $E
## [1] 0.5
adj_mat = rbind(c(0,3,1,0,0,0,0,0,0), #A's neighbors
c(3,0,2,2,0,0,0,0,0), #B's neighbors
c(1,2,0,0,2,0,0,0,0), #C's neighbors
c(0,2,0,0,1,0,1,1,0), #D's neighbors
c(0,0,2,1,0,2,0,2,0), #E's neighbors
c(0,0,0,0,2,0,0,0,0), #F's neighbors
c(0,0,0,1,0,0,0,1,0), #G's neighbors
c(0,0,0,1,2,0,1,0,1), #H's neighbors
c(0,0,0,0,0,0,0,1,0) #I's neighbors
)
rownames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
colnames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
# Convert adjacency matrices to igrpah objects for all three graphs.
ig = graph.adjacency(adj_mat, mode="undirected", weighted=TRUE, add.colnames = "name")
print(ig)
## IGRAPH 5ee3d84 UNW- 9 12 --
## + attr: name (v/c), weight (e/n)
## + edges from 5ee3d84 (vertex names):
## [1] A--B A--C B--C B--D C--E D--E D--G D--H E--F E--H G--H H--I
adjacency_matrix = list(adj_mat)
perms = mle.getPermMovie_memory(subset.nodes = c("A", "B"), ig, output_filepath = "/Users/lillian.rosa/Downloads/CTD/vignette/biased_walker", movie=TRUE, zoomIn = FALSE)
## [1] "Calculating permutation 1 of 2..."
## [1] "Calculating permutation 2 of 2..."
# Get perms as list object, with no images generated
perms = mle.getPerms_memory(S = c("A", "B"), G)
## [1] "Calculating permutation 1 of 2."
## [1] "Calculating permutation 2 of 2."
# This network walker is less sensitive, but node permutations can be calculated ahead of time (along the lines of a dynamic programming solution),
# so thousands of subsets can be reasoned about very quickly.
# Generate PNGs to animate the network walker.
perms = mle.getPermMovie_memoryless(subset.nodes = c("A", "B"), ig, output_filepath = "/Users/lillian.rosa/Downloads/CTD/vignette/unbiased_walker", movie=TRUE, zoomIn = FALSE)
## [1] "Calculating permutation 1 of 2..."
## [1] "Calculating permutation 2 of 2..."
# Get perms as list object, with no images generated
S = c("A", "B")
perms = list()
for (n in 1:length(S)) {
perms[[n]] = mle.getPermN_memoryless(n, G, S, misses.thresh = log2(length(G)))
}
## [1] "Calculating permutation 1 of 9."
## [1] "Calculating permutation 2 of 9."
names(perms) = S
We’re going to go back to our data using the Arginase deficiency network model, and the Miller et al (2015) dataset. ## IV.I Choose your node subset.
# Maximum subset size to inspect
kmx=15
# Some global variables from before.
adjacency_matrix = list(arg_icov)
G=vector(mode="list", length=length(V(ig_arg)$name))
G[1:length(G)] = 0
names(G) = V(ig_arg)$name
# Get our node subset associated with the $KMX highest perturbed (up or down) in our first Arginase deficiency sample.
S_arg = sort(abs(arg_data[,1]), decreasing=TRUE)[1:kmx]
print(S_arg)
## phenylcarnitine homoarginine n-acetylarginine
## 7.549122 6.717854 6.224659
## x - 12339 x - 12681 4-guanidinobutanoate
## 6.064197 6.007678 5.448070
## 3-methyl-2-oxovalerate leucine 4-methyl-2-oxopentanoate
## 5.390644 5.116677 4.517897
## ornithine x - 11381 isoleucine
## 4.457695 4.197347 4.075734
## hippurate x - 19330 arginine
## 3.848954 3.725342 3.406121
# Get the adaptive network walker's paths starting from each node in the subset S_arg.
perms = mle.getPerms_memory(S = names(S_arg), G)
## [1] "Calculating permutation 1 of 15."
## [1] "Calculating permutation 2 of 15."
## [1] "Calculating permutation 3 of 15."
## [1] "Calculating permutation 4 of 15."
## [1] "Calculating permutation 5 of 15."
## [1] "Calculating permutation 6 of 15."
## [1] "Calculating permutation 7 of 15."
## [1] "Calculating permutation 8 of 15."
## [1] "Calculating permutation 9 of 15."
## [1] "Calculating permutation 10 of 15."
## [1] "Calculating permutation 11 of 15."
## [1] "Calculating permutation 12 of 15."
## [1] "Calculating permutation 13 of 15."
## [1] "Calculating permutation 14 of 15."
## [1] "Calculating permutation 15 of 15."
ptBSbyK = mle.getPtBSbyK_memory(names(S_arg), perms)
ptID = colnames(arg_data)[1]
data_mx.pvals=apply(arg_data[,which(colnames(arg_data) %in% diagnoses$id)], c(1,2), function(i) 2*pnorm(abs(i), lower.tail=FALSE))
res = mle.getEncodingLength(ptBSbyK, t(data_mx.pvals), ptID, G)
2^-res[,"IS.alt"]
## [1] 1.615509e-03 8.077544e-04 4.038772e-04 2.019386e-04 1.009693e-04
## [6] 5.048465e-05 2.524233e-05 7.888227e-07 7.703347e-10 1.925837e-10
## [11] 9.629183e-11 4.814592e-11 2.407296e-11 1.203648e-11 6.018240e-12
2^-res[,"d.score"]
## [1] 1.000000e+00 6.474570e-03 6.292531e-05 8.176575e-07 1.329077e-08
## [6] 2.597783e-10 5.934544e-12 2.481370e-12 3.743036e-11 2.454003e-12
## [11] 8.864218e-14 3.498948e-15 1.498843e-16 6.929281e-18 3.435771e-19
data_mx = arg_data[,which(colnames(arg_data) %in% diagnoses$id)]
data_mx = data_mx[,1:8]
ptBSbyK = list()
for (pt in 1:ncol(data_mx)) {
ptID=colnames(data_mx)[pt]
S_arg = sort(abs(data_mx[,pt]), decreasing=TRUE)[1:kmx]
perms = mle.getPerms_memory(S = names(S_arg), G)
ptBSbyK[[ptID]] = mle.getPtBSbyK_memory(names(S_arg), perms)
}
res = list()
t = list(ncd=matrix(NA, nrow=ncol(data_mx), ncol=ncol(data_mx)),
jacdir=matrix(NA, nrow=ncol(data_mx), ncol=ncol(data_mx)))
rownames(t$ncd) = colnames(data_mx)
colnames(t$ncd) = colnames(data_mx)
rownames(t$jacdir) = colnames(data_mx)
colnames(t$jacdir) = colnames(data_mx)
for (i in 1:kmx) {
res[[i]] = t
}
for (pt in 1:(ncol(data_mx)-1)) {
print(sprintf("Patient %d vs...", pt))
ptID=colnames(data_mx)[pt]
for (pt2 in (pt+1):ncol(data_mx)) {
print(sprintf("Patient %d.", pt2))
ptID2=colnames(data_mx)[pt2]
tmp = mle.getPtSim(ptBSbyK[[ptID]], ptID, ptBSbyK[[ptID2]], ptID2, data_mx, NULL)
# Similarity based on mututal information, by subset size.
# plot(1:kmx, tmp$NCD)
# Similarity based on Jaccard set similarity with directionality, by subset size.
# plot(1:kmx, tmp$dirSim)
for (k in 1:kmx) {
res[[k]]$ncd[ptID, ptID2] = tmp$NCD[k]
res[[k]]$jacdir[ptID, ptID2] = tmp$dirSim[k]
res[[k]]$ncd[ptID2, ptID] = tmp$NCD[k]
res[[k]]$jacdir[ptID2, ptID] = tmp$dirSim[k]
}
}
}
for (pt in 1:ncol(data_mx)) {
for (k in 1:kmx) {
res[[k]]$ncd[pt, pt] = 0
res[[k]]$jacdir[pt, pt] = 0
}
}
res_all = res
# Multi-dimensional scaling
# if you have diagnostic labels associated with the colnames(data_mx), send them using diagnoses parameter
res_ncd = lapply(res_all, function(i) i$ncd)
res_jacdir = lapply(res_all, function(i) i$jacdir)
ncd = mle.getMinPtDistance(res_ncd)
jacdir = mle.getMinPtDistance(res_jacdir)
diags = colnames(data_mx)
diags[which(diags %in% diagnoses$id[which(diagnoses$diagnosis=="Argininemia")])] = "ARG"
diags[which(diags %in% diagnoses$id[which(diagnoses$diagnosis!="Argininemia")])] = "negCntl"
names(diags) = colnames(res_all[[1]]$ncd)
# alpha weighs infromation from NCD with information from Jaccard-dir. alpha must be between 0 and 1.
p = htmltools::tagList()
i = 1
for (alpha in seq(0, 1, 0.1)) {
patientSim=alpha*ncd+(1-alpha)*jacdir
colnames(patientSim) = colnames(res_all[[1]]$ncd)
rownames(patientSim) = colnames(res_all[[1]]$ncd)
p[[i]] = as_widget(plot.mdsSim(patientSim, diags, k=2, NULL))
i = i + 1
# Hierarchical clustering
# A PNG image called "ptSimilarity.png" will save in the path set as parameter "path" in plot.hmSim function call.
plot.hmSim(patientSim, path=getwd(), diags)
# K-NN
plot.knnSim(patientSim, diags, diag="ARG")
}
## [1] "Tp = 4, Tn= 1, Fp = 3, Fn=0"
## [1] "Sens = 1.00, Spec= 0.25"
## [1] "Tp = 4, Tn= 1, Fp = 3, Fn=0"
## [1] "Sens = 1.00, Spec= 0.25"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
## [1] "Tp = 4, Tn= 2, Fp = 2, Fn=0"
## [1] "Sens = 1.00, Spec= 0.50"
p
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels